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The two-stream instability has been mooted as an explanation for a range of astrophysical appli- 
cations from GRBs and pulsar glitches to cosmology. Using the first nonlinear numerical simulations 
of relativistic multi-species hydrodynamics we show that the onset and initial growth of the instabil- 
ity is very well described by linear perturbation theory. In the later stages the linear and nonlinear 
description match only qualitatively, and the instability does not saturate even in the nonlinear case 
by purely ideal hydrodynamic effects. 
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I. 



INTRODUCTION 



The two-stream instability occurs generically when two coupled, inter-penetrating fluids have a sufficiently large 
relative velocity. The expectation is that the instability arises when a perturbation appears to move in different 
directions with respect to each fluid. Originally studied in magnetized plasmas modeled by the collisionless Vlasov- 
Maxwell equations [U [2], where a detailed understanding of the full mode spectrum exists (see e.g. [3 ), it is only 
recently (see [4 ) that it has been studied in the context of pure relativistic hydrodynamics. The well-understood 
plasma instability has astrophysical applications in GRBs (see e.g. [3 ) and the pulsar emission mechanism ([HE]), 
whilst the hydrodynamic instability has been suggested as a key mechanism behind pulsar glitches ([71|8]) and is of 
relevance in the cosmological context [9 . 

Most existing explorations of the two-stream instability have focused on its development in the linear regime, as - 
in astrophysically interesting applications such as pulsar glitches - the instability is expected to saturate at relatively 
low amplitude. However, the mechanism by which the instability saturates is not clear. The obvious possibility is 
that non-ideal effects such as shear viscosity will stop the growth. However, given the difficulties in constructing a 
consistent and stable relativistic non-ideal theory this avenue can only be pursued phenomenologically. An alternative 
possibility is that the instability saturates due to nonlinear effects within the framework of the ideal theory, possibly 
via the generation of internal shocks converting the unstable modes to heat. A final possibility that is unlikely to be 
generically applicable is that the fluids may be dynamically driven out of the instability window, perhaps by external 
forces or non-local effects, as happens in the cosmological case (j9]). 

In this paper we will study the nonlinear development of the relativistic hydrodynamic two-stream instability in 
simplified cases using numerical simulations. For reasons detailed later we are not able to consider shock propagation. 
However, the simulations can investigate if shock formation is a possibility, or whether the high frequency oscillations 
dominate, at which point we would expect non-ideal effects to be important. Our results suggest that when the 
instability is triggered by small perturbations the two-stream instability grows until the solution is dominated by high 
frequency oscillations. By comparing with linearized time-domain solutions we see that the nonlinear coupling has 
only a small effect, and is insufficient to saturate the instability. 



We consider the system of relativistic multiple fluids as introduced by Carter in [10 and detailed in the review of 
Andersson and Comer [11]. Here the notation largely follows [11], assumes the existence of a spacetime metric gab of 
signature — h ++, uses Roman letters from the start of the alphabet - a^b^c^ . . . - as 4-spacetime indices, and from 
the middle - i^j^k^ . . . - as 3-space indices in the 3 + 1 split. The characters X, Y, Z will be used as labels indicating 
the different fluids, or different species, which will not be implicitly summed over except where explicitly stated. Units 
where the speed of light c = 1 are used throughout. 



The basic fluid quantities are the number density currents where X is a species label. Here we do not consider 
reactions or particle creation, meaning that the number densities are conserved. Hence the currents obey the continuity 
equations 



The system is closed using the master function A. The equations of motion - the Euler equations - follow by 
minimizing the action defined using this master function. The master function is defined in terms of all possible scalar 
invariants n^y = —Qah^yJ^Y' shorthand notation n\y^ = n^- is often used. The Euler equations are written in 
terms of the conjugate momenta /i^, defined by 



II. THE SYSTEM 



A. 



General form 



Van^ = 0. 



(1) 




(2) 



Using these definitions the Euler equations follow as 



= 2n^V[af^f] = 0. 



(3) 
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B. 3 + 1 decomposition 

As a test-bed we will consider flat spacetimes in standard plane symmetric Cartesian coordinates. Using (t, i) to 
represent the time and spatial coordinates, the continuity equation ([T]) becomes 

dtn'^ + difi], = (4a) 

and the Euler equation (|3| can be written as 

j 

dtiif - diiif - 2^%4 = 0. (4b) 



C. Numerical implementation 

The nonlinear numerical simulations solve the equations of motion Q. Although these describe conservative ideal 
(multi-species) hydrodynamics, they are not written in balance law form. Therefore discontinuous "solutions" will 
not necessarily be the correct (entropy satisfying) solution, irrespective of the numerical techniques employed - the 
derivation of a suitable entropy satisfying balance law form is required future work before multifluid systems with 
shocks can be studied. For simplicity we will therefore use a basic finite difference discretization of the equations 
of motion Q and note that the numerical solutions should not be trusted in two regimes: firstly, when the spatial 
gradients are sufficiently large to appear discontinuous on the grid, and secondly when the solution is dominated by 
high frequency (relative to the numerical grid) components. 

The implementation used here relies on the Method of Lines where the time integration uses either the standard 
third order strict stability preserving (SSP) Runge-Kutta method, or the standard fourth order non-SSP Runge-Kutta 
method. The spatial discretization uses central differencing of either second or fourth order. We have tested whether 
the addition of Kreiss-Oliger dissipation (of either third or fifth order respectively) makes any difference, and it does 
not. In all simulations we set the timestep by imposing a CFL constant At/ Ax of 0.25, and have checked that the 
results are insensitive to the precise value. 

The numerical solution of Eq. ^ gives solutions for and /i^. This is sufficient to describe the system completely, 
but in order to numerically compute the terms for the next update it is necessary to construct the spatial components 
of the number density ffux, n^. To do this we note that the definition of the conjugate momenta, Eq. ([2|, can be 
seen as a system of nonlinear algebraic equations for the n^. Solving this nonlinear algebraic system gives the correct 
value for the scalars n^, from which we can compute the master function A and its derivatives. This allows Eq. ([2| 
to be viewed as a linear system for the spatial components n^, as required. 

III. LINEARIZED SOLUTION 

To study the time domain behaviour of the multifiuid system, and as a comparison test for the nonlinear simulations, 
we construct linearized solutions of two coupled fiuids on a 2L periodic domain [— L, L] in one spatial dimension. The 
restrictions on the number of fiuids and spatial dimensions can be dropped in the following analysis at the cost of 
substantially complicating the explicit calculation of the final solution. 

A. Linearized system 

Start by restricting the equations of motion Q to one spatial dimension x, which after linearizing gives the system 

dtSn^x + dxSn^ = 0, (5a) 
dtS^^ - dj^f = 0. (5b) 

Note that, from the definition of the conjugate momenta given in Eq. ([2|, we have (see pT] ) 

5^^=M^^Snl (6) 

where here, and in the future, Z is an abstract species index to be summed over, and where the matrices can be given 
explicitly (e.g. [11^ section 11.3]). 
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We note that the linearized system can be written as 

dtSnt^ + djn^ = 0, 



(7a) 
(7b) 



We will not use this form in the linearized solution, but will use Eq. ^ to construct numerical solutions, using the 
techniques of Sec. |II C| for direct comparison with the nonlinear results. 



B. Transformed system 

We now use a Fourier-Laplace analysis, firstly taking the discrete Fourier transforms by writing 

fcmax 



^^x = ^x,[/c]exp[ia;/ex] 

/c=0 

where cj/^ = iik/L. Then the Laplace transform 



leads to the fully transformed equations 



IWfe 



(8) 

(9) 

(10a) 
(10b) 



Solution 



Combining the results from Eq. (10), the solution for follows from the linear system 



where the matrix A is given by 



and the vector V by 



X,[fe] 



tt ^'^Z,[m]- 



(11) 

(12) 
(13) 



The solution for A^^ follows from Eq. (10a). Given these solutions we can invert the Laplace and Fourier transforms 
to construct Sn^ at any future time. 



D. The instability 

We note that the qualitative change from a stable two-fluid system to an unstable one can easily be seen with this 
linearized solution. In the formal solution given by Eq. ([iTJ the inverse of the matrix A^^ can be written 



(A 



xz\ 



"114=1 (s - i'^fe^'i) 



(14) 



where the denominator, a quartic polynomial in s, is the determinant of A^^. This form is chosen so that the roots 
Vi are independent of the frequency ujk- On inverting the Laplace transform to construct the solution for Sn^ we 
find that the time dependent behaviour is encoded in exponentials of the form exp(ia;/er^t). Therefore, the resulting 
linearized solution is stable only if the roots have vanishing imaginary part and, when unstable, the growth rate 
is linearly related to the frequency as expected. This precisely mirrors previous calculations (e.g. [4J) where the 
stability was usually found more straightforwardly in terms of the dispersion relation. 
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Entrainment 


1 


1 ^ ^ 
^22 


8 
5 


9 
5 


N/A N/A 




Chemical coupling 





2-^/2 1 


1 


4 
3 


p 1.1 1.1 






TABLE I. The parameters used in the simulations below, based on the master function given in Eq. (17). In each case one 
parameter, indicated with p, controls the coupling between the two species. In the pure entrainment case where K12 = the 
values of ax are irrelevant. 



Species 1 at rest, entrainment case, k^=0.5 



Species 1 at rest, cliennical coupling case, k^^=Q.025 

0.2 1 ^ ^ ^ ^ 




0.2 0.4 0.6 0.8 1 
Bacl^ground velocity species 2 



0.2 0.4 0.6 0.8 1 
Background velocity species 2 



FIG. 1. The two-stream instability acts when the roots of the determinant of a particular matrix (see Eq. ( 14 ) and accompanying 
text) have non- vanishing imaginary part. The left plot shows the entrainment case, and the right the chemical coupling case 
- the full parameters are outlined in the text and Table ^ For this particular choice there is little qualitative difference in the 
instability window. 



IV. RESULTS 



In what follows we shall always use an equation of state inspired by [9^ and JV2! and encompassing the key expected 
behaviour. First note ([12j) that the Lorentz factor of one fluid in the frame of the other is 



"12 

nin2 



(1 



implying that A, defined by 



n2 



encodes the velocity difference between the two fluids. Then we use a master function with general form 

2 



A(n?,n^,n?2) 



x=i 



(15) 



(16) 



(17) 



We would approximately expect the /^a term to encode the entrainment effects, and the hzi2 term to encode the 
chemical coupling. 

We look at the entrainment and chemical coupling cases separately, insofar as this separation makes sense at the 
nonlinear level. The precise parameters employed are given in Table ^ In the pure entrainment case the strength of 
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FIG. 2. Convergence tests of the nonlinear code against the hnearized solution. In both cases the nonlinear code is fully fourth 
order accurate and the curves are scaled to show the fourth order convergence. The entrainment case is shown, and in both 
cases the perturbation is a sine wave period 27r, amplitude S = 10~^, in m. In the left panel the fluids are uncoupled, — 0, 
and both are at rest in the background. In the right panel the fluids are coupled with At a =0.5 and the first fluid also has a 
non-trivial background velocity of 0.1 whilst the second is at rest. 



the coupling is determined by the parameter a^a- For the pure chemical coupling case we consider parameters relevant 
to cosmology (see [9 ), where the strength of the coupling is determined by the parameter ni2. The chemical coupling 
case also checks that the methods work in the limit mx ^ as expected. Fig. [T] shows the instability window for 
these representative cases where the background is given by the simple choice n\ = 1. 

The initial data chosen here is always written as a perturbation about a background, even in the nonlinear case. 
We write (in 1 + 1 dimensions) 

= nxWx (-1, ^x, 0, 0) , ^x = (l - 4) (18) 

where Wx is the Lorentz factor for the X species. Perturbations about the background to either nx or vx (or both) 
have been tested. In what follows we shall concentrate on perturbations in nx, as no qualitative difference between 
the cases has been found. In the majority of the cases we focus on sinusoidal perturbations of amplitude 5^ with 
different frequencies compatible with the size of the domain. These show clean behaviour at very low numerical 
resolution and are easy to study. Other perturbations (such as the Gaussian shown below) have been tested with no 
qualitative change in the key features, but may be less flexible (because of the periodic boundaries) and may require 
higher resolution (due to the formation of steep gradients without forming discontinuities). 



A. Convergence 



Firstly we use the linearized solution outlined in section III to benchmark the nonlinear code in the stable regime. 



This is only possible for a certain range of initial perturbation amplitudes, outside of which the convergence of the 
error is limited by nonlinear effects or floating point precision. In Fig. [2] we show simple convergence tests for a 
single mode perturbation in the entrainment case - similar results are seen in the chemical coupling case. The master 
function parameters are as detailed above. In addition, for the uncoupled case we have tv/\ = ^ and for the coupled 
case we choose tz^ = 1/2. The uncoupled case is the simplest possible - both fluids are at rest in the background 
with nx = 1. In the coupled case a velocity difference is imposed by giving the first fluid a velocity of 0.1. The 
perturbation is imposed only in the first fluid and is a simple sine perturbation, period 27r, amplitude 5 = 10~^, in 




FIG. 3. A direct comparison of the time evolution of the number density components in the coupled case used in the convergence 
plot in Fig. [2] The linearized solution is given by the solid black lines, and the nonlinear numerical evolution (using only 32 
points) is given by the blue circles. With these parameters we see that the fluids are fully coupled within one crossing time, 
but that the profiles retain the general shape given by the initial perturbation. 



the number density ni. Similar results are seen when the initial perturbation is a Gaussian, although considerably 
higher resolution is required to resolve the spatial gradients. 

We note that the results converge as expected for low resolutions. At higher resolution we typically see results that 
are not perfectly convergent, either due to the tiny nonlinear couplings starting to dominate over the linear effects, 
or because the numerical error is affected by floating point precision. Where the dominant effect is the nonlinear 
couplings we still see the expected 5 e(f- convergence of the nonlinear code. 

A direct comparison of the components of the number densities is given in Fig. [3j Over the period of approximately 
one crossing time we see the excellent agreement between the linearized solution and the nonlinear evolution, even 
at extremely low numerical resolution. The fluids rapidly couple with this choice of entrainment parameter - within 
this time, which is approximately one crossing time, the amplitude of the second fluid is comparable to that of the 
first which contained the perturbation. The form of the perturbation is also retained through the evolution. 



B. Instability growth 



We first consider the entrainment case with = 1/2. As an illustration we look at a case where the first fluid is 
perturbed with a simple sine wave as before, but is at rest, whilst the second fluid is unperturbed but has a velocity 
of 0.6 in the background so that the two- stream instability acts. Representative results are shown in Figs. [4] and [5j 
We note first that the space-time development of the instability appears to be (a) dominated by the high frequency 
components, (b) not spatially localized to the perturbation (as shown in particular when the initial perturbation is a 
narrow Gaussian, as in Fig. [5|, and (c) closely following the linearized solution, where we assume that the linearized 
solution is "seeded" by constant frequency low amplitude perturbations in all modes in addition to the perturbation 
used in the numerical simulation. 
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, -6 unstable case -^ t = , -6 unstable case -^ t = 0.25 
X 10 xJO xMO xJO 




-1 -0.5 0.5 1 -1 -0.5 0.5 1 -1 -0.5 0.5 1 -1 -0.5 0.5 1 



FIG. 4. A direct comparison of the time evolution of the number density components in the coupled case within the instability 
window for the two-stream instability. The linearized solution is given by the solid black lines, and the nonlinear numerical 
evolution (using 128 points - every fourth shown) is given by the blue circles. This is the entrainment case with At a = 0.5 and 
a background velocity difference of 0.6. The instability grows rapidly in both linear and nonlinear case; to compare, the linear 
solution seeds all modes up to /c = 64 with small random noise. In the early stages the behaviour closely follows the stable case 
shown in Fig.js] By t — 0.5 (bottom left panel) the instability is visible in the high-frequency oscillations for the linear solution 
(solid line), but is difficult to see in the nonlinear numerical solution (circles). By i — 0.75 (bottom right panel) the instability 
is clear in both the nonlinear numerical solution and the linear solution, where the full extent of the oscillations no longer fit 
within the scale of the plot. The two-stream instability grows fastest in the highest frequency modes, seeded by hand in the 
linear case and by numerical error in the nonlinear case. The exponential growth is not spatially localized and the coupling 
shows the instability in all components of all species. 



1. Time- frequency behaviour: linear case 

More detail can be seen when the solution is studied in time-frequency space. These results, however, mix numerical 
and nonlinear effects. In order to disentangle these effects as far as possible we first consider solely the linearized case, 
as shown in Fig. [6j In the linear case we expect the different frequency modes to behave independently. The change 
in power of a given frequency mode with time will therefore illustrate only effects due to the two-stream instability 
where relevant and the numerical method employed. The linearized solution, as constructed in Sec. [Ill C[ is shown 
in the left panel and shows simple characteristic behaviour, with power-law behaviour in the instability growth time 
scale. This behaviour is not replicated in the numerical simulations of the linear system, shown in the right panel. 
This must be purely due to the numerics used. This numerical effect can be modelled analytically. 

Following the work of Lele [13] we note that the effect of the finite difference scheme, when approximating spatial 
derivatives, is equivalent to replacing the exact Fourier transform relation dxf iujf with dxf ^ i^'/, where uj' 
encodes how accurately the numerical method captures modes of a given frequency. For the centred fourth order 
differencing that we focus on here a straightforward calculation (or appropriately applying the general results of [13]) 
gives 

c^'(^) = (19) 

Here we use uj = Trcj/Zcmax to scale to Lele's units. Repeating the linear analysis using uj^ in place of u throughout leads 
to the results in the central panel of Fig. [6| The match between these adjusted "exact" results and the results from 
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FIG. 5. A direct comparison of the time evolution of the number density components in the coupled case within the instability 
window for the two-stream instability. The linearized solution is given by the dotted black lines, and the nonlinear numerical 
evolution (using 512 points - every point shown) is given by the blue dashed line. This is the entrainment case with k,a = 0.5 
and a background velocity difference of 0.6, as in Fig. [4] but here the initial perturbation is a Gaussian, hence the higher 
resolution required. The instability grows rapidly in both linear and nonlinear case; to compare, the linear solution seeds all 
modes up to k — 256. 



numerical simulation in the right panel of Fig. [6] is extremely good. The remaining minor differences are likely due 
to the choice of using uniform amplitude "noise" to seed the analytic calculation at higher frequencies, as compared 
to the numerical simulation where the noise likely appears from floating point round-off error. This assumption is 
supported by the amplitude of the seed noise used in this comparison, which is ^ 10~^^ x S: the relative floating point 
accuracy of the linearized perturbation. 



2. Time-frequency behaviour: nonlinear case 



We next sketch our expectations of how the nonlinear effects would modify the behaviour. Considering a single fluid 
and taking solely the discrete Fourier Transform of the nonlinear system of Eq. Q, we would expect the nonlinear 
modes to satisfy coupled equations of the qualitative form 

Stnf,]+A[,](*)nf,] =0, (20) 

where ® is the (cyclic) convolution product. In the linear case the matrix A^i^^ would only contain the constant (zero 
frequency, k = 0) term. In the nonlinear case the terms depend on the master function and the data - in other words, 
on the non-trivial h^j^y 

With initial data dominated by say two frequencies ki and k2 the convolution couples modes at the harmonics with 
frequencies aki + bk2 with a, b integers. For the typical initial data used here the dominant frequencies will be from 
the background {k = 0) and the initial perturbation with frequency /cq, leading to harmonics at integer multiples of 
k(). This effect can be clearly seen in numerical nonlinear simulations as shown in Fig.[7[ In this case the background 
is chosen so that the velocity difference is 0.1 so the two-stream instability does not act. The fluids couple and the 
higher harmonics are excited, but only the frequencies associated with the background, the initial perturbation and 
the first harmonic would be noticeable in the spatial snapshots. 
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FIG. 6. A time-frequency plot of the growth of the instabiUty in the hnearized case. Both time and frequency are scaled to be 
resolution independent - the numerical simulation shown uses 1024 points (/cmax = 512). The norm of the complex component 
of the discrete Fourier Transform of 6n\ is shown - the result for other components and quantities is qualitatively the same. 
The left plot shows the behaviour for the linearized solution. As expected the highest frequencies blow up the fastest. The 
central plot shows the behaviour for the linearized solution, adjusted as in Eq. (19) to take into account the behaviour of the 
spatial differencing scheme. This would be the expected behaviour from a numerical evolution of the linearized equations. The 
right plot shows the results from the numerical linear evolution. We see that the adjustment for the spatial differencing scheme 
completely explains the numerical behaviour, suppressing the growth of the instability at high frequencies. 



Once the background is modified to allow the two-stream instability to act we can see both the numerical suppression 
of the high frequency modes (as shown in the linear case in Fig.[6| and the nonlinear coupling of the higher harmonics 
(as shown in the stable case in Fig. [7|. Fig. [s] repeats the calculation exciting higher harmonics as in Fig. [oj but the 
background is modified so that the velocity difference is 0.6 and the two-stream instability acts. Whilst all modes 
are excited to some degree by numerical error, only those that are coupled via harmonic overtones are excited to 
a level that is numerically significant (i.e. greater than 10~^^ here). The figure then shows precisely these excited 
modes growing exponentially as expected, initially independently of the other modes, with the suppression at higher 
frequency qualitatively matching the behaviour seen in the linear case. 

We can now look at the original case, studied for the linear problem in Fig. [Oj Here the only mode initially excited 
is the lowest mode on the grid. As seen in Fig.|9j there is still a qualitative match between the adjusted linear solution 
and the numerical evolution, even in the nonlinear case. The "bulk" features, including the approximate growth rates 
and the suppression of the growth at high frequencies, remain the same. Whilst there are also noticeable differences 
at late times, there is no clear pattern in either the space-time development or the time-frequency growth, except for 
the increased coupling between the different frequencies at late time. 



C. The instability onset 



The parameters chosen for both the entrainment and chemical coupling case have a "window" where the two-stream 
instability acts, as shown by Fig. [l] By modifying the relative velocity of the background in the initial data, we can 
check that the onset of the instability occurs at the same point with the nonlinear code. 

In both the entrainment and chemical coupling case we find that, when using small perturbations of amplitude 
S = 10~^ as above, the onset of the instability (with increasing relative velocity) is, to numerical precision, identical 
to that predicted from the linear analysis. In the chemical coupling case we can also check the upper edge of the 
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Numerical power in 6n^^: k^^^ = 64, 128, 256. Initial frequency %l^2 




Numerical power in 6n^^: k^^^ = 64, 128, 256. Initial frequency 3:r/64 




FIG. 7. A time- frequency plot showing the harmonic couphng expected in the nonUnear case. As the grid resolution (and 
hence frequency maximum /cmax) is varied the initial perturbation is modified to give a constant frequency with respect to the 
grid. In all cases the perturbation has a single frequency. The nonlinear coupling immediately excites all the higher harmonics, 
although only the lowest harmonic (at twice the frequency of the initial perturbation) is excited at a level that would be visible. 
The initial parameters are chosen such that the fluids couple but the two-stream instability is not acting. 
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Numerical power in 6n^^: k^^^ = 64, 128, 256. Initial frequency %IZ2 




max max max 



Numerical power in 6n^^: k^^^ = 64, 128, 256. Initial frequency 37r/64 




max max max 



FIG. 8. A time- frequency plot of the growth of the instabiUty where the initial frequency (relative to the grid) of the perturbation 
is held fixed. As the grid size is increased the same number of harmonics are excited by the nonlinear coupling, as in the stable 
case shown in Fig. [71 Each excited harmonic then grows exponentially due to the two-stream instability, with the pattern 
resembling (at early times) the linearized case shown in Fig. [g]). As in the stable nonlinear case in Fig. [t] there is no visible 
nonlinear coupling except to the immediate harmonics. In the top row the overtones are always at higher frequencies. The 
choice of period in the bottom row means that some overtones at lower frequencies are excited, meaning that in the bottom 
left plot all frequencies on the grid are excited by nonlinear couplings. 
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FIG. 9. A time-frequency plot of the growth of the instabihty in the nonhnear case. This figure is the nonlinear equivalent 
of Fig. [6] The linearized and adjusted solutions are seeded with noise from the first timestep of the numerical solution; this 
is noticeably larger than in the linear case. The qualitative comparison between the adjusted and numerical solution is still 
reasonable, explaining the suppression of the growth at high frequencies, but the nonlinear coupling does modify the solution 
significantly. 



window; that is, the nonlinear results show the instability for 0.29185 < < 0.69985, and are stable otherwise. In 
the entrainment case we are unable to check the upper edge of the window (which is at Av ^ 0.958) as the velocities 
required are too large for the nonlinear code to successfully evolve. 

Finally we attempted to produce initial data that would generate shocks before the two-stream instability becomes 
important, to see whether the nonlinear couplings could ever be expected to be more important than the instability 
growth. In our experiments, only extremely large initial perturbations {S ~ 10~^), combined with a choice of master 
function parameters such that the two- stream instability growth is as slow as possible, show signs of characteristic 
breaking before the two-stream instability sets in. Even in these cases it is possible for the two-stream instability to 
dominate simply by increasing the grid size, and so admitting grid frequencies that grow sufficiently rapidly. Within 
this purely ideal hydrodynamics case using the shearing box approximation it does not seem possible to suppress the 
instability. 



D. Higher dimensions 

All the results so far have been restricted to 1 + 1 dimensions in the shearing box (periodic boundaries) approxi- 
mation. We have performed numerical simulations that retain the shearing box approximation in higher dimensions 
- in particular, detailed comparisons in 2 + 1 dimensions were calculated. No qualitative differences were found. This 
can be understood by transforming to the frame of the background flow for one species. In this frame, there are only 
two possible effects from the additional spatial dimension. 

First, the angle between the perturbed flow and the background velocity difference will change the coupling. This 
is indeed the case: our simulations confirm that the projection of the perturbation onto the velocity difference must 
be non-zero for the instability to act. However, in this (generic) case, the behaviour of the instability growth is 
qualitatively identical to the 1 + 1 dimensional case. 

Second, the coupling between the species should potentially lead to a change in the alignment between the species, 
modifying the relative velocity between them. This could potentially change the growth of the instability. However, 
explicit simulations performed by modifying the angle of the perturbation with respect to the relative flow showed 
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no visible difference in, for example, the growth of A^. Thus it does not seem possible for purely local hydrodynamic 
effects, even at the nonlinear level, to modify the behaviour of the instability. 



V. DISCUSSION 



The two-stream instability has been mooted as an explanation for a range of astrophysical applications from GRBs 
and pulsar glitches to cosmology. We have used numerical simulations to study the nonlinear development of this 
instability when the species are modelled using coupled relativistic hydrodynamics. 

Our simulations show that the onset of the instability when only the local, purely hydrodynamic behaviour is 
considered, perfectly matches the predictions of linear theory. The restricted analysis in 1 + 1 dimensions is sufficient, 
provided the initial data is appropriately interpreted, when there is no change in either the spacetime, the external 
forces, or the coupling parameters. When these restrictions are relaxed, for example in the cosmological case considered 
by [9 where the expansion of the universe is modelled by varying the coupling parameters, then the instability growth 
can be curtailed. 

The initial growth of the instability is also well described by linear theory, after adjustment to account for the 
known properties of the numerical simulations. At the very late stages where the instability dominates is there a 
noticeable difference between the linear and nonlinear behaviour, which has no obvious pattern. 

This work is a necessary first step towards generic nonlinear simulations of relativistic multifluids. As multifluid 
effects such as entrainment are expected to be important in, for example, heat conduction [TT, charge conduction and 
resistivity [15 , and neutron superfluids [16 , high accuracy nonlinear simulations incorporating detailed microphysics 
must include multifluids. The main technical barrier to incorporating these effects in current simulations is the lack 
of either a balance law form or a non-conservative form that explicitly controls entropy at discontinuities. Our future 
work is aimed at overcoming this technical hurdle before combining the multifluid effects with nonlinear relativistic 
elasticity [17 and appropriate numerical techniques for interfaces [18 to produce such simulations. 
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